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We present an analysis of the Rayleigh-Plesset equation for a three dimensional vacuous bub¬ 
ble in water. In the simplest case when the effects of surface tension are neglected, the known 
parametric solutions for the radius and time evolution of the bubble in terms of a hypergeometric 
function are briefly reviewed. By including the surface tension, we show the connection between 
the Rayleigh-Plesset equation and Abel’s equation, and obtain the parametric rational Weierstrass 
periodic solutions following the Abel route. In the same Abel approach, we also provide a discussion 
of the nonintegrable case of nonzero viscosity for which we perform a numerical integration. 
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INTRODUCTION 

It is well established that the size evolution of unstable, spherical cavitation bubbles is governed by the Rayleigh- 
Plesset (RP) equation [1—3] 

P w (rR + =p- Poo - + 2p w R S j . (1) 

In 0 p w is the density of the water, R(t) is the radius of the bubble, p and P ^ are respectively the pressures 
inside the bubble and at large distance, a is the surface tension of the bubble, and p w is the dynamic viscosity of 
water. In the simpler form with only the pressure difference in the right hand side, equation m was first derived by 
Rayleigh in 1917 [l| but it was only in 1949 that Plesset developed the full form of the equation and applied it to the 
problem of traveling cavitation bubbles Q • In the second half of the last century a steady progress has been achieved 
with driving forces from engineering, medical, sonoluminescence, microfluidics, and even pharmaceutical applications 
of the cavitation phenomena Q. In addition, the interest in the analytical and numerical solutions of cavitation 
dynamics remained considerable as these solutions can lead to a better control and understanding of the bubble 
collapse processes. In an effort to discern the peculiar features of the usual three-dimensional collapse, Prosperetti [1] 
and more recently Klotz [&] worked out generalizations to TV-dimensional bubble dynamics. 

In this paper, we first review the analytic solutions of the RP equation in terms of hypergeometric functions when 
the surface tension is neglected, and in terms of Weierstrass elliptic functions when the surface tension is taken into 
account jj]. In the latter case, we employ an Abel equation approach which is a novel mathematical way of looking 
to the nonlinear evolution of cavitation bubbles. On the other hand, when the viscous term is introduced, we show 
that an Abel equation with a non constant invariant occurs. Since it is not yet known how to find analytical solutions 
to this equation (if any), we resort to numerical integration of the stiff RP equation from t = 0 to t = t c , where t c is 
the time of collapse of the bubble. 


SIZE EVOLUTION AND COLLAPSE OF A SPHERICAL BUBBLE WITHOUT SURFACE TENSION: 

REVIEW OF CANONICAL RESULTS 

We first consider the idealized case whereby the viscosity of the water is neglected, since p w -C 1, and we further 
discard the effect of surface tension a. For this case, there are also recent analytical approximations in Q, which are 
further discussed in 9f, but here we are concerned with the standard results. 

Let us consider a vacuous p = 0 bubble of radius R which is surrounded by an infinite uniform incompressible fluid, 
such as water, that is at rest at infinity. We remark that ‘infinity’ in the present context refers to distance far enough 
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away from the initial position of the bubble, and we further assume that the pressure at infinity is constant, = 
const. Neglecting the body forces acting on the bubble, we have from equation © 


2 RR + 3 R 2 = -2 


. P„ 


( 2 ) 


Since R 2 R is an integrating factor of © consequently we obtain by one quadrature 

r 3 r 2 = j ? 3 + e . 

3 Pw 

Using the initial conditions i?(0) = Ro and i?(0) = 0 we find the integration constant to be 


p _ 2 Poo „3 

6 “ 3—^0 ’ 
Pw 


and hence, we obtain 


R 2 


2 P^ 

3 Pw \RJ 


(3) 


Note that one can find a simple novel particular solution for R(t) by substituting © into © to obtain the Emden- 
Fowler equation 


R = At n R m 


(4) 


with A = — n = 0, m = —4, and particular solution 


25 P^R 3 / >/6 [fn[ R \ 5 
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Rp(t) = 
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t + 
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However, this solution is obtained under the assumption of a nonzero integrating factor Q and therefore it fails 
to satisfy the second initial condition. Other solutions that do not satisfy the initial conditions have been found 
previously by Amore and Fernandez [§]. As noticed in [l()j |, equation © can be also viewed as a conservation law for 
the dynamics of the radius of the bubble, since its kinetic energy can be expressed as 


2 7T Pw R 3 R 2 = (i?g - R 3 ) . (6) 

To proceed with the integration of equation © we will use the set of transformations as given by Kudryashov 
namely R = S e ,dt = R s dr, where e, 5 are constants that depend on the dimension of the bubble, and S,t are the 
new dependent and independent variables, respectively. Applying the transformations upon ©, we obtain the new 
dynamics in S and r 


S 2 = l^^(R 3 S~ 3e - 1)S 2+2e5 ~ 2e . (7) 

3 Pw ^ 

To find S, one sets e = A, anc } 8 = N + 1, where N = 3 is the dimension of the bubble, which will in turn reduce © 
to the simpler equation 


S T = 



By integrating the above with S'(O) = Rq we obtain the rational solution 


5(r) 


R 3 o 

®r 2 + 1 ’ 
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(9) 
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where for convenience we set 23 = ^-Rq = § Once we determine S, we can find the parametric solutions for 

the bubble radius R(t) and evolution time of the bubble t(r) 16] 


R(t) = 


R o 


(23r 2 + l)i 


t(r)=R 4 0 [ 
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dt; 
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o (®C 2 + 1) 3 


The integral for the evolution of the time for bubble can be calculated analytically in terms of hypergeometric functions 
to give |6| 


Hr) = ^ 
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To achieve the time of collapse one needs to allow r 00 . This leads to lim,--^ i(r) = 0.000908681 if we use, in 
S.I. units, \p w \ = 1000 kg/m 3 , [i?o] = 10 -2 to and [Poo] = 101325 Pa. We point out that the same asymptotic value 
for the collapse time has been also obtained, albeit using a different approach, by Obreschkow et al. [8[ and it is also 
obtained from Eq. [IT] below by direct integration. Once we solve for r as a function of R from the first equation of 
m, and substituting it into the second equation of m we obtain the closed-form solution 


t(R) — Rq 
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1 4 3 
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which is plotted as R{t) in Fig. [T] 
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FIG. 1: Radius of the bubble in the absence of surface tension according to equation (1121) for \p w ] = 1000 kg/m 3 , [Ho] = 10~ 2 m 
and [Poo] = 101325 Pa. 



Next, we will find the time for the total collapse t c of the bubble by integration of equation (]3|), which yields 

rH o ^> 3/2 


t c = 


Rn 


1 /3 pu 
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This is just the N = 3 case of the general IV-dimensional formula given in [j|. If we let R = Rosin 2 ^ 3 d, where 
6 £ [0, 7t/ 2] then the integral (fill) transforms to 


tc 


2R 0 /37^ r/ 2 
3 \2P x J a 


sin 2/3 9 dO, 




















4 


which by comparing with the integral relation for Beta function, namely, 


B(m, n) = 2 


/*7r/2 


COS 


sin 2 ” -1 9 dO 


leads to 


t c 


Ro /3a7 r (1 5\ 
3 V 2 P x \2’ 6j 




0.914681 Ra 



(14) 


This solution can be also obtained from (J4]) by using the parametric solution 2.3.1 —2.2 from Polyanin’s book 11]. If we 
insert in (1T41) the same S.I. numbers as used in the hypergeometric asymptotics, we obtain again t c = 0.000908681 sec. 


SURFACE TENSION INCLUDED VIA ABEL’S EQUATION 


When the surface tension term added, equation m can be written in the form 


3 

2 r: 




where we define K\ = 


-P 00 --P 


and K 2 = —. 

Pw Z Pw 


Proceeding as in [12], we first show that solutions to a general second order ODE of type 

R + f 2 (R)R + h(R) + fi(R)R 2 + fo(R)R 3 = 0 
may be obtained via the solutions to Abel’s equation 1171) of the first kind (and vice-versa) 

dy 


dR 


= MR) + h{R)y + h(R)y 2 + h{R)y 3 


using the substitution 

R = r)(R(t)) , 

which turns (1161) into the Abel equation of the second kind in canonical form 

vv + fs{R) + f 2 (R)y + h(R)y 2 + fo(R)y 3 = 0 . 

Moreover, via the inverse transformation 

rim)) = 1 


(15) 


(16) 


(17) 


(18) 


(19) 


( 20 ) 


y{R(t)) 

of the dependent variable, equation (1T51) becomes (1171) . The invariant of Abel’s equation, see [l3|, can be written as 


m) = 


1 (df 2 


fs ~ f2~r^ - /1/2/3 d 77/2 


3 \dR dR ' g- 

tt Abel’s equation is integrable. 

By identification of equation © with ©, we see that /i(l?) = f 2 {R) = f 0 (R) = 0, and f 3 {R) = 


( 21 ) 


and when is a constant is an indication that Abel’s equation is integrable. 

JL — n f~( m — K 1 _i_ Aa 

R w R 2 , 

and hence the Kamke invariant is $(1?) = 0, therefore Abel’s equation (1171) becomes the Bernoulli equation 


j^ = h(R)y + h(R)y 3 


( 22 ) 


which, by one quadrature, has the solution 

y{R) = 


±V?> Ri 


V3D - 2KJF - 3K 2 R 2 ' 


(23) 











By using equations (l20l) and (fl 8 l) . we obtain 
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R z = 


32) - 2K 1 R 3 - 3K 2 R 2 
3 R 3 


and via the same initial conditions we obtain the integration constant 

V = ^-(3K 2 +2K 1 R 0 ) , 


(24) 


which gives 


R 2 


2K l (R 3 - R 3 ) + 3K 2 (R 2 - R 2 ) 
3 R 3 


(25) 


Notice that when a = 0 —> I\ 2 = 0, and p = 0 —> AT — ^, then the above becomes equation ([3]). The new energy 
with surface tension is 


2 tt Pw r 3 r 2 = — [Pco(i?o - R 3 ) + MRl - R 2 )] ■ 

Thus, at the time of collapse we have Rq(PooR 0 + 3 a) = 0.42443 J for a surface tension [a] = 10 -3 N/m. 

To integrate equation (1251) . first let us write it in a more convenient way, as 


(26) 


• 2 03 Cll 

R ~& + R +a ° 


(27) 


with coefficients defined as 03 = (K 2 + ^j x Rq'), a\ = —K 2 , and ao = and we will use the same set of 

transformations, namely R = S e , dt = R s dr which give in turn 


Si = 


s 2+2e6 


(aoS~ 2e + ai S- 3E + a 3 S~ be ) 


(28) 


where now we set e = —3/N = —1, and 6 = ^-±1 = 2. Thus, we obtain the Weierstrass elliptic equation 


which in standard form is 


via the linear substitution 14] 


Sz = ao + aiS + a 3 S 3 


Pr =4p 3 - g 2 p - g 3 


S'(r) = — p(r + 29.709; g 2 , g 3 ) . 
03 


(29) 


(30) 


(31) 


The germs of the Weierstrass function g 2 ,g 3 are given by 

a 3 a 3 _ K 2 Rf 
4 ~ 4 

a 0 a§ _ K, R,q 
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) Z - 


16 24 v 3 

Substituting AT and K 2 into the germs, the solution to equation (l29l) becomes 


S(T) = WpBR^)' v{t + 29 ™; *■») ■ < 33 ) 

where g 2 = 3.37751 • 10 ” 11 m 8 /sec 4 , g 3 = 1.92645 • 10 -8 m 12 /sec 6 , and the constant in the front of Weierstrass 
function from (l33l) takes the value of 59215.2 sec 2 /m 5 . Once S is known we can find the parametric solutions for the 
radius and time of the bubble with surface tension as 
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R(t) = 
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S{t) 6 p w 

T di _ RtiPooRo + 3a) 2 
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Related plots are presented in Fig. [2] 
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FIG. 2: Parametric solutions for time of evolution (top), radius of the bubble (middle), and radius vs time (bottom) from 
equation (1341) when surface tension is present for [pw] = 1000 kg/m 3 , [Ro] = 10~ 2 m, [-Poo] = 101325 Pa and [<r] = 10~ 3 N/m. 
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FIG. 3: Comparison between the exact solutions without (dotted curve) and with surface tension (dashed curve). Notice that 
when we have surface tension the time of collapse is smaller. 


THE RAYLEIGH-PLESSET EQUATION WITH VISCOSITY 


When we include the viscosity, equation CD reads 

R + f 2 (R)R + f 1 (R)R 2 + f 3 (R) = 0, (35) 

where f 2 (R) = and K 3 = is a constant, and fi(R), f 3 (R) being the same functions as before. Thus, Abel’s 
equation 03 becomes 

^ = f\(R)y + f 2 (R)y 2 + h(R)y 3 • (36) 

Now, we also find the Karnke invariant according to equation m and we obtain 


m) = 


K 3 [4A'| - 9R{3K 2 + hK\R)\ 


In terms of the physical variables of the system, the invariant is 

2y w [6Ay 2 w - 9p w R(6a + 5PooR)} 


*{R) = 


27p3i?6 


(37) 


(38) 


and because is not a constant, we will try to reduce (15C>1) using the Appell invariant instead. 

First, we will eliminate the linear term via the transformation y(R) = R? z(R) to obtain the reduced Abel equation 


dz 

dR 


h 2 {R)z 2 + h 3 (R)z 3 , 


(39) 


where h 2 (R) = =, and h 3 (R) = {K\R + K 2 )R. 

According to the book of Karnke [l3[, for equations of the type (15^1) for which there is no constant invariant one 
should change the variables according to 


which lead to the canonical form 


where 


z(R) = z(C(R)) , 

C(#) = [ h 2 (R)dR = 2K 3 VR , 


dz 

dC 


vf(C)5 3 


*(C) 


h 3 (R{ 0) 
h*(R( 0) 


C 3 (b 3 + b 5 ( 2 ) 


(40) 


(41) 


(42) 








is the Appell invariant, and the constants bi are 


K-i = vpj 
8 K* 2 4 > 4 ’ 

Ki _ PqqpI, 

32 K* ~ 2^a ■ 


(43) 


Choosing a dynamic viscosity of water \fi w \ = 1.002 cP = 1.002 • 10 3 kg/(m ■ sec ) the bi constants take the values 
of 63 = 9.68789 • 10 14 sec 2 /m 5 , while 65 = 7.63836 • 10 32 sec 4 /m 10 . The units of £ are to® /sec. This Abel equation 
is not integrable through quadratures, but numerically we integrate the RP equation (l35l) . from t = 0 to the point of 
stiffness which is the point in time where the bubble collapses, see the numerical solution on Fig. |4] 


R(t) 



FIG. 4: Numerical solution for the RP equation (1351) when surface tension and viscosity are both present for [p,„] = 1000 kg/m 3 , 
[Po] = 10~ 2 m, [Poo] = 101325 Pa, [a] = 10 -3 N/m and [g w ] = 1.002 • 10 ~ 3 kg/(m ■ sec). 



FIG. 5: Comparison between the numerical solution (continuous curve) of Fig. [4] and the exact solutions of Fig. [3]. As one 
can see, the bubbles collapse fastest when there is surface tension (dashed curve) followed by the no surface tension and no 
viscosity case (dotted curve). The slowest time is achieved when we have viscosity together with surface tension (continuous 
curve). Right panel is the zoomed in version of the figure in the left panel. 


CONCLUSION 

In this work, we have considered the RP equation for the size evolution of a bubble in water. In the first part of 
the paper, we have surveyed the standard results in the absence of surface tension but in a different way from those 
usually pursued in the literature. We have obtained the closed form hypergeometric solutions of Kudryashov and 
Sinelshc.hikov although in a different but equivalent form. From an Emden-Fowler form of the RP equation in this case 
we have also obtained a particular solution which however does not satisfy the second initial condition. In the presence 
of surface tension and viscosity, we have employed a new approach based on Abel’s equation. When only the surface 
tension is included, we have obtained the known parametric rational Weierstrass solutions, whereas when viscosity 
is added, the corresponding Abel equation does not have a constant invariant, which explains the nonintegrability in 
this case. A numerical integration obtained from this nonintegrable Abel route is presented graphically. 
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